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ABSTRACT 

This paper describes the calculation of flutter stability 
characteristics for a transonic forward swept fan configuration 
using a viscous aeroelastic analysis program. Unsteady Navier- 
Stokes equations are solved on a dynamically deforming, body 
fitted, grid to obtain the aeroelastic characteristics using the 
energy exchange method. The non-zero inter-blade phase 
angle is modeled using phase-lagged boundary conditions. 
Results obtained show good correlation with measurements. It 
is found that the location of shock and variation of shock 
strength strongly influenced stability. Also, outboard stations 
primarily contributed to stability characteristics. Results 
demonstrate that changes in blade shape impact the calculated 
aerodynamic damping, indicating importance of using accurate 
blade operating shape under centrifugal and steady 
aerodynamic loading for flutter prediction. It was found that 
the calculated aerodynamic damping was relatively insensitive 
to variation in natural frequency. 

INTRODUCTION 

Reducing noise has become a significant criterion in the 
design of modem turbomachinery fans. For fans with 
supersonic tip speeds, shock waves play a significant role in the 
generation of noise. Aggressive blade designs with forward 
sweep and lean are being used to improve the noise 
characteristics. These geometric features coupled with the 
presence of shock wave can lead to flutter instability. Flutter 
problems, usually detected during the design phase, result in 
program delays and cost overruns, adversely impacting the 
development and maintenance costs. Numerical flutter 
prediction methods will help in achieving the design objectives 
and developing a flutter free design. 

Several methods have been developed for the prediction of 
aeroelastic characteristics of turbomachines with various 
degrees of fidelity. Aeroelastic analyses based on energy 
exchange between vibrating blades and the surrounding fluid, 
have been reported for turbomachines using semi-analytical 
methods (Lane and Friedman [1], Smith [2]), linearized 
methods (Hall and Clark [3], Verdon [4]), Euler methods (He 
[5], Gerolymos and Vallet [6], Bakhle, et al. [7]) and viscous 
methods (Giles and Haimes [8], Siden [9], He and Denton 
[10]). A limited number of coupled aeroelastic analyses of 
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turbomachine configurations have also been reported [11-13]. 
Williams, Cho, and Dalton [11] used a linear panel method to 
solve the eigenvalue problem. Gerolymos [12] and Srivastava 
and Reddy [13] solved the coupled aeroelastic equations based 
on an inviscid aerodynamic analysis. The majority of the 
methods are either analytical, applicable to two-dimensional 
geometry and/or ignore the viscous effects. 

The modem fans being designed to reduce noise are swept 
forward as well as leaned in the plane of rotation. These 
features coupled with transonic flow result in aeroelastic 
problems. For turbomachinery blade rows, flutter is observed 
primarily in single natural mode. Single mode flutter can be 
captured using the energy exchange method, which calculates 
the flow of energy between the vibrating blade and the 
surrounding fluid. A three-dimensional viscous aeroelastic 
analysis code TURBO-AE, based on energy exchange method 
has been developed and reported in [14]. The objective of the 
present study is to investigate the stability of a forward swept 
transonic fan geometry. The transonic fan analyzed, showed 
flutter at part speed. The TURBO-AE code is applied to 
calculate the flutter observed in the wind tunnel and to 
understand the flow features that caused the instability. 

Modem fan blades, by virtue of having a small thickness 
ratio coupled with forward sweep could be highly flexible in 
the tip region. High rotational speed and large aerodynamic 
loading result in significant changes in blade shape. It is 
important to accurately account for these changes in the 
analysis. To investigate the sensitivity of damping calculations 
to changes in steady blade deflection, calculations are also 
carried out by deforming the blade at the tip to change the blade 
twist by 0.5 deg such that the incidence angle at the section is 
increased. Also, results are obtained to investigate the impact of 
time step used in the analysis and sensitivity to natural 
frequency variation. These studies are carried out to gain an 
understanding of issues that need to be addressed in flutter 
analyses of similar fan configurations. 

THE TURBO-AE CODE 

The aeroelastic solver TURBO-AE is briefly described in 
this section. The solver can model multiple blade rows 
undergoing harmonic oscillations with arbitrary inter-blade 
phase angles (IBPAs). It is based on a Navier-Stokes unsteady 
aerodynamic solver for internal flow calculations of axial flow 
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turbomachinery components, TURBO [15,16]. Viscous effects 
are modeled using the Reynolds-averaged Navier-Stokes 
equations. The two-equation k-£ turbulence model is used for 
closure. The aerodynamic equations are solved using a finite 
volume scheme. Flux vector splitting is used to evaluate the 
flux Jacobians on the left-hand side. The right-hand side fluxes 
are discretized using the high order Total Variation 
Diminishing (TVD) scheme based on Roe's flux difference 
splitting. Newton subiterations are used at each time-step to 
maintain accuracy. Symmetric Gauss-Seidel iterations are 
applied to the discretized equations for improved convergence. 

The aeroelastic characteristics of the rotor are obtained by 
calculating the energy exchange between the vibrating blade 
and its surrounding fluid. Positive work on blade indicates 
instability. The aeroelastic analysis is performed by first 
obtaining the “steady” aerodynamic solution for a given 
operating condition. The blades are then forced into a 
prescribed harmonic motion (specified mode, frequency, and 
IBPA) to calculate the unsteady aerodynamic response and 
thereby work-per-cycle W, Eq. 1. 

— > 

W=( J) l-p.dA»Adt (1) 

surface 

Where, p is the surface pressure, A is the area of the blade 

surface and X is the blade displacement associated with the 
forced harmonic vibration. The work W is then used to 
calculate the aerodynamic damping y, using the average kinetic 
energy K E of the blade over one vibration cycle [17], Eq. 2. 


The blade motion is simulated using a dynamic grid 
deformation technique. Phase-lagged boundary conditions are 
used to calculate the non-zero IBPA vibrations, which 
eliminates the need to model multiple blade passages. 

RESULTS AND DISCUSSION 

A scale model of an experimental transonic fan with 22 
blades (Fig. 1) was tested in a rig. The design operating 
condition for the model tested was mass flow of 98.89 lbs/s 
with a relative tip Mach number of 1 .4. The fan fluttered at part 
speed in the first natural mode [18], between the operating line 
and predicted stall line. The first natural mode is shown in 
Fig. 2, hotter colors signifying larger displacements. The mode 
shows a high degree of twist bend coupling, especially in the 
outboard section.. 

The results are presented here for 90 percent rotational 
speed. A grid with approximately 213,000 grid points per blade 
passage, 46 axial by 47 radial points on each surface, and a 
uniform tip-gap is used in the analysis. The blade shape and 
modal properties for 100 percent speed are used. 

Figure 3 shows the steady aerodynamic loading on the 
blade surfaces for 90 percent speed at a condition near the 
operating line. Shock wave structures are clearly evident from 
these figures. An oblique shock appears near the leading edge 


on the suction surface over the entire span. A normal shock 
extends across the blade passage in the aft section of the blade 



Figure 1. Transonic fan 



on the suction surface. For the in-board sections, this shock is 
outside the blade passage (in front of blade leading edge) on the 
pressure surface side. However, the shock is completely 
ingested within the passage in the outboard sections of the 
blade. This shock pattern is typical for a flow condition 
referred to as moderate to high loading for blades with a 
supersonic leading edge operating in subsonic axial flow. As 
the back pressure is increased, moving the operation point 
towards the stall line, the shock moves forward and becomes 
stronger. This results in flow separation on the suction surface 
in the middle section of the blade. The separation is large 
enough to stall the blade resulting in a break down of the 
numerical analysis. The observed breakdown of the analysis 
due to flow separation could not be verified, as velocity profiles 
were not measured during the experiment. 
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Figure 3. Steady blade surface pressure near operating line for 
90 percent speed 

Dynamic Analysis 

Experiments showed blade flutter in the first natural mode, 
Fig. 2, (natural frequency 351 Hz) for 32.73 deg IBPA (2 Nodal 
Diameter forward traveling wave) [18]. To establish the least 
stable mode the analysis was carried out for several pressure 
ratios and lBPAs. 

Flow separation for higher back pressures was observed in 
the hub region of the blade resulting in unsteadiness of the 
mean flow. Because of the mean flow unsteadiness, one needs 
to ensure that the blade vibration amplitude is large enough so 
as to minimize the impact of mean flow unsteadiness on the 
unsteadiness due to blade vibrations, yet should be small 
enough to maintain linearity. The calculated time history of 
mass flow is shown in Fig. 4 for flow with and without blade 
vibration for a back pressure of 16.2 pounds per square inch 
(psi). The blade is forced into vibration at 180 deg IBP A and 
351 Hz. It can be seen from Fig. 4 that for the vibration 
amplitude chosen for the analysis, the unsteadiness due to blade 
vibration is at least an order of magnitude larger than the mean 
flow unsteadiness. It can also be seen that the flow converges to 
periodicity in roughly 10 oscillations. 



Figure 4. Variation of mass flow for the transonic fan, with and 
without blade vibration 


Figure 5 shows the aerodynamic damping calculated for 
several different back pressures and IBPAs. For lower back 
pressures, the least stable IBPA was found to be 0 deg. 
However, as the back pressure was raised moving the fan 
operating condition towards the stall line, the 32.73 deg IBPA 
(2 nodal diameter forward traveling wave) became the least 
stable. Flutter was observed in the wind tunnel for 32.73 deg 
IBPA. 



Figure 5. Variation of aerodynamic damping with back 
pressure and IBPA for the transonic fan 


The calculated variation of aerodynamic damping with 
back pressure is shown in Fig. 6 at 32.73 deg IBPA. As may be 
seen in the figure, as the back pressure is increased, the 
aerodynamic damping decreases rapidly, dropping to 
approximately 0.2 percent of critical damping at a back 
pressure of 16.4 psi. Increasing the back pressure above 16.4 
psi for 90 percent speed resulted in stalled flow with large 
separation on the suction surface in the mid section of the 
blade. This separation was found to be shock induced and 
prevented a steady mean flow from which the blade vibration 
analysis could be carried out. These results indicate that the 
32.73 deg IBPA is the least stable IBPA for first natural mode, 
as observed in experiments, and that the aerodynamic damping 
rapidly decreases with increasing back pressure. 

From Fig. 6 it is evident that the characteristics of variation 
of aerodynamic damping changes significantly as the back 
pressure is varied. Initially the damping is not affected much 
by an increase in back pressure, at least up to the operating line, 
which falls close to a back pressure of 15.7 psi. Further 
increase of the back pressure results in a dramatic change in 
damping variation, becoming inversely proportional to 
variation in back pressure. To understand the reason for such a 
change, flow details were investigated for back pressure values 
of 15.2, 15.9 and 16.4 psi. Figure 7 shows the distribution of 
work on the blade surface. Majority of the blade shows little 
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Figure 6. Variation of aerodynamic damping for 16.4 PSI back 
pressure. Forced vibration in the first natural mode at 351 Hz 
and 32 deg 1BPA 


blade positive work is being done on the blade indicating a 
destabilizing contribution. These plots, when superimposed 
with relative Mach number plots, indicate that the positive 
work regions are centered around the shock location. 

To investigate this in more detail, the pressure coefficient 
at 90 percent span, around the shock location on the suction 
surface, is plotted in Fig. 8 for the three back pressures. 
Instantaneous pressure coefficient at five instants of time over 
one period of oscillation is shown. The steady pressure 
coefficient is also shown to highlight the variation of shock 
strength from the mean position. It is clear that as the back 
pressure is increased the shock moves upstream. Also, the 
variation in the pressure change associated with the shock, 
changes significantly. It is these changes that contribute to 
positive work on the blade. These figures indicate that the 
location of shock and associated changes in strength are 
responsible for the change in characteristics of the variation of 
damping and are responsible for flutter. 


contribution towards establishing the flutter characteristics of 
the blade. Outboard sections of the blade contribute most 
significantly. This has a side benefit from an analysis 
perspective. Resolving the flowfield on the inboard sections is 
probably not significant and hence coarse grid may be used in 
this region, concentrating more grids in the outboard sections of 
the blade. These figures also show that over certain areas on the 
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Figure 7. Work distribution on blade surface for various back 
pressures 
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Figure 8. Pressure coefficient variation over one period of 
oscillation cycle for back pressures of 15.2, 15.9, 16.4 psi 


The blade natural vibration frequency for a rotor is 
proportional to the rotational speed. Hence, for 90 percent 
speed the blade natural frequency would be lower than that at 
100 percent speed. To investigate the impact of change in 
frequency, a sensitivity study was carried out by varying the 
natural frequency. Calculations were performed for 20 percent 
higher as well as 20 percent lower natural frequency for 100 
percent speed at 32.73 deg IB PA. These changes in frequency 
resulted in changes of less than 0.1 percent of critical damping 
ratio, indicating a low sensitivity to natural frequency. It 
should be noted here that changes in rotational speed will also 
impact the mode shape, however, in the present study no 
attempt was made to alter the mode shape. 

Analysis was also performed at the second natural 
vibration mode, which was stable in the experiment. The 
aerodynamic damping calculated was much higher than those 
for the first mode, indicating the second mode to be more stable 
than the first mode. These results show that the analysis 
successfully calculates the natural mode and IBPA of instability 
identifying the flutter characteristics of the transonic fan. 

For the back pressure of 16.4 psi, the experiment showed 
flutter, whereas the analysis predicted the fan to be marginally 
stable. The analysis correctly predicted the least stable mode, 
however, it did not predict a negative aerodynamic damping. A 
possible reason for that was thought to be the use of an 
inaccurate operating blade shape. The analysis at the 90 percent 
speed was performed using the design speed geometry and 
characteristics. To investigate the influence of changes in blade 
shape on the aerodynamic damping, the analysis was carried 
out by altering the twist distribution, such that the blade twist 
changed from mid-span to tip varying linearly from zero 
change at mid span to 0.5 deg change at the tip. The change in 
twist distribution is such that it increases the incidence of the 


blade. Structural analysis showed the blade twist change at the 
tip to be of this order for the 90 percent speed. The variation in 
aerodynamic damping with oscillation cycles is shown in 
Fig. 9. The changes in twist reduce the damping moving it 
closer to instability. Earlier studies have shown smaller time 
steps to have moderate impact on the damping calculation, see 
Fig. 6. As the calculated damping is close to zero, the impact of 
smaller time steps was also investigated. Most of the analysis 
was carried out using 100 steps per oscillation cycle in order to 
minimize the computational cost. The computational cost 
varies linearly with the number of steps per cycle used in the 
analysis. Three different time steps were used, 100, 200 and 
300 steps per oscillation cycle. The impact of time step size is 
also shown in Fig. 9. From 100 steps per cycle to 200 steps per 
cycle there is a measurable change in damping. Whereas using 
300 steps per cycle the change is small and does not justify the 
additional computational cost. From this study it was decided 
that the 100 steps per cycle is good for investigating the trends 
of the behavior. For cases where a more accurate damping 
value is required, 200 steps per cycle should be used. Since the 
difference is approximately 0.2 percent of critical damping, use 
of smaller time steps will be significant only in the 
neighborhood of zero damping or flutter boundary. 

Figure 10 compares the steady pressure coefficient and 
work contributions from the suction and pressure surfaces at 
the 95 percent span location. The results for the original grid 
are (Fig. 10(a)) obtained using 100 steps per cycle, where as for 
the deformed grid the results (Fig. 1 0(b)) are for 200 steps per 
cycle. These comparisons are to understand the differences in 
calculated damping for the deformed shape. Primarily the 



Figure 9. Effect of time steps and change in twist distribution 
for a 0.5 deg change in twist at tip of the blade 
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suction surface contributes to the positive work whereas the 
pressure surface contributes to negative work providing a 
stabilizing influence. The area of positive work is centered 
around the normal shock on the suction surface. However, there 
are two main differences that are observed from these figures. 
First, we find that the shock on the suction surface has moved 
upstream by roughly 5 percent. In addition, the work 



(a) Original Geometry 


contribution associated with this shock has also moved. Only 
minor changes in the level of work are observed. A second 
change, which is more significant from the work perspective, is 
seen on the pressure surface. The area of negative work in the 
leading edge region is reduced, and the peak is higher but 
narrower for the deformed shape. This change contributes 
mainly to a reduction in the damping for the deformed shape. 
The primary reason for this difference is the motion of the 
shock on the pressure surface. Figure 11 shows the steady 
pressure contours and Fig. 12 show the instantaneous pressure 
contours at five instants of time over the vibration cycle for 
steady and deformed geometry at 95 percent span. The steady 
contours show that a normal shock wave extends from the 
suction surface across the blade passage. For the original grid, 
the shock intersects the pressure surface at the leading edge, 
whereas for the deformed geometry the shock is just ahead of 
the leading edge. Figure 12 shows the shock movement over 
the oscillation cycle for the two geometries considered here. 
Instantaneous pressure contours near the pressure surface 
leading edge are shown for five time instants over one cycle of 
vibration. For ease of comparison, these contours are plotted 
next to each other going from left to right with increasing time. 
The pressure surface, leading edge, and time instants 
corresponding to the contour are marked on the figures. From 
Fig. 12 a it can be seen that the shock stays on the blade surface 
for the original geometry, whereas for the deformed geometry 
(Fig. 12 b), the shock intersects the blade for part of the 
vibration cycle, moving on and off the blade surface over the 
oscillation cycle. This shock motion on and off the blade 
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(b) Deformed geometry 


Original geometry 


Deformed geometry 


Figure 10. Pressure coefficient and work distribution at 
95 percent span 


Figure 11. Pressure contours for 16.4 psi at 95 percent span 
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surface is the reason for differences in work contribution on the 
pressure surface between the two geometries. The effect of this 
motion is to reduce the overall stabilizing effect from the 
pressure surface, reducing the aerodynamic damping for the 
deformed geometry. 




Figure 12. Instantaneous pressure contours over one vibration 
cycle with time period T, near pressure surface leading edge at 
95 percent span for back pressure of 16.4 psi 

Approximately 95 Mega Words of computer memory and 
one hour of computational time per 100 steps was required on 
the Cray C-90 computer located at the NAS facility of NASA. 


CONCLUDING REMARKS 

An aeroelastic analysis program based on the Navier- 
Stokes equation has been used to calculate the flutter 
characteristics and to investigate the behavior of aerodynamic 
damping variation for a forward swept transonic fan geometry. 
Calculated variations of aerodynamic damping with back 
pressure and IB PA is presented for a condition where flutter 
occurred during testing. Although a negative aerodynamic 
damping was not calculated, the analysis yielded some good 
results and trends. The analysis correctly predicted the mode 
and I BP A of flutter. Also, the trend for calculated aerodynamic 
damping clearly indicated a negative aerodynamic damping 
would result with further increasing the back pressure. 
Unfortunately, the analysis showed that a stalled flow field 
would emerge for any further increase of the back pressure. 
Although flow separation was captured along with vortex 
shedding, deep blade stall posed a numerical problem. Results 
also showed that a majority of the blade surface does not 
contribute towards the stability characteristics. But the 
stabilizing and destabilizing effects were concentrated near the 
shock location and were restricted to outboard stations of the 
blade. 

Detailed examination of the flow features indicated that 
normal shock on suction surface has a destabilizing effect and 
results in a large area of positive work centered around the 
normal shock. The area of the destabilizing effects depended 
on the location of the shock. Higher back pressure moved the 
normal shock forward increasing the area of positive work. 
This coupled with the mode shape resulted in an increased 
variation of the blade surface pressure over one oscillation 
cycle causing the destabilizing influence to increase. 

A study was also performed to understand the effect of 
various parameters on the calculation of aerodynamic damping. 
The study showed that the aerodynamic damping was 
insensitive to the variations in vibration frequency. However, it 
should be noted that only the frequency was altered without 
altering the mode shape. It was also found that the blade shape 
under operating conditions is important. A small change of 0.5 
deg in twist at the blade tip led to a better correlation of 
damping with measurements, where the original shape had 
predicted marginal positive damping. Changes in blade shape, 
although moved the suction surface shock further upstream in 
the outboard sections, do not measurably change the level of 
positive work. On the other hand, the negative work on 
pressure surface was reduced in the front section of the blade 
reducing the overall blade stability, resulting in better 
correlation with measurements. Thus it is important to use the 
correct shape in the analysis for more accurate calculation of 
damping, especially near the flutter boundary. 

Changes in time step for unsteady analysis were found to 
have a moderate impact. Although the trends did not change, 
the smaller time step reduced the aerodynamic damping. 
However, these differences would be of significance for 
damping calculations only in the vicinity of neutral stability. 
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From these studies the following conclusions can be drawn 

that could be useful in similar analyses: 

(a) Shock location and its strength have a strong influence 
on blade stability 

(b) In-board sections of the blade do not contribute 
towards establishing the stability characteristics, thus a 
coarser grid in this region can be used to reduce 
computational cost. 

(c) Accurate blade operating shape is important for 
accurate damping prediction as it establishes the shock 
structure on the blade, which in turn affects the 
stability. Variations in vibration frequency were found 
to be less significant. 

(d) Larger time steps can be used away from the neutral 
stability region for determining the trends, however, 
close to neutral stability; smaller time steps need to be 
used for accurate stability predictions. 
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